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Matlab Software for Iterative Methods and 
Algorithms to Solve a Linear System 


PROF. D.A. GISMALLA 


Abstract — The term "iterative method" refers to a wide 
range of techniques which use successive approximations to 
obtain more accurate solutions .In this paper an attempt to 
solve systems of linear equations of the form AX=b, where A is a 
known square and positive definite matrix. We dealt with two 
iterative methods namely stationary (Jacobi, Gauss-Seidel, SOR) 
and non-stationary (Conjugate Gradient, Preconditioned 
Conjugate Gradient).To achieve the required solutions more 
quickly we shall demonstrate algorithms for each of these 
methods .Then using Matlab language these algorithms are 
transformed and then used as iterative methods for solving these 
linear systems of linear equations. Moreover we compare the 
results and outputs of the various methods of solutions of a 
numerical example. The result of this paper the using of 
non-stationary methods is more accurate than stationary 
methods. These methods are recommended for similar situations 
which are arise in many important settings such as finite 
differences, finite element methods for solving partial 
differential equations and structural and circuit analysis . 


Index Terms —Matlab language, iterative method, Jacobi. 


I. Stationary iterative methods 

The stationary methods we deal with are Jacobi iteration 
method, Gauss-Seidel iteration method and SOR iteration 
method. 

A. Jacobi iteration method 

The Jacobi method is a method in linear algebra for 
determining the solutions of square systems of linear 
equations. It is one of the stationary iterative methods where 
the number of iterations is equal to the number of variables. 
Usually the Jacobi method is based on solving for every 
variable Xi of the vector of variables X - =(x 1 ,x 2 ,....,x n ) , 
locally with respect to the other variables .One iteration of the 
method corresponds to solve every variable once .The 
resulting method is easy to understand and implement, but the 
convergence with respect to the iteration parameter k is slow. 

B. Description of Jacobi’s method: - 

Consider a square system of n linear equations in n variables 

AX=b (1.1) 

where the coefficients matrix known A is 

A =(«i,;) ij - 1(1) n 

The column matrix of unknown variables to be determined X 
isA' : = (r L ,i a . 

and the column matrix of known constants b is 


= Gq, b 2 f ... f b n ) 

The system of linear equations can be rewritten 

(D+R)X=b (1.2) 

where A=D+R , D =(&,- J e = l(l)?a 

is the diagonal matrix D of A and R = L + U , 
where L and IT are strictly lower and 
strictly Upper matrix of A 

Therefore, if the inverse D -1 exists and Eqn.(1.2) can be 
written as 

X=D -1 (b-RX) (1.3) 

The Jacobi method is an iterative technique based on solving 
the left hand side of this expression for X using a previous 
value for X on the right hand side .Hence , Eqn.(1.3) can be 
rewritten in the following iterative form after k iterations as 
Xt^sZT 1 (b-RJC**- 1 )) ,k=l,2,.... (1.4) 

Rewriting Eqn.(1.4) in matrix form . We get by equating 
corresponding entries on both sides 

(iL 1 iff n (k— 0 

x i ~ . 2 - 1 — = l q !J x j 

LL “■ j*L 

i=l,...,n ,k=l,2,.... (1.5) 

We observe that Eqn. (1.4) can be rewritten in the form 

= TX&-V + C k=l,2,., (1.6) 

where T=D _1 (-L-U) and C = D ~~b 

or equivalently as in the matrix form of the Jacobi iterative 

method 

X^=D- 1 C— L— UjX^+D-'b k=l,2,.... 
where T=D -1 (-L-U) and C = D”'b 

In general the stopping criterion of an iterative method is 
to iterate until 




< E 


for some prescribed tolerance E> 0 .For this purpose, any 
convenient norm can be used and usually the norm , i.e. 

|| X -k>_ x ^-0|| 

— <e,X k * 0. 


This means that if X^ 1 is an approximation to X ^'Jhen 
the absolute error is |x ^ - X ^ _1J [ , and the relative error is 

—- provided that X ■ ^ 0. This implies that a 


vector X 1 can approximate X to t significant digits if t is the 

||X-X'IL 

largest non negative integer for which _L p 7 j—- <5*10 “■ 

C. The Matlab Program for JACOBI 

The Matlab Program for Jacobi’s Method with it’s 
Command Window is shown in Fig.(1.1) 
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D. Gauss-Seidel iteration method 

The Gauss-Seidel method is like the Jacobi method, 
except that it uses updated values as soon as they are available 
.In general, if the Jacobi method converges, the Gauss-Seidel 
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method will converge faster than the Jacobi method, though 
still relatively slowly. 

Now if We consider again the system in Eqn.(l.l) 

AX=b (1.1) 

and if We decompose A into a lower triangular component 
L* and a strictly upper triangular component U. Moreover if 
D is the diagonal component of A and L is the strictly lower 
component of A then 

A= L*+U 

Where 


L* —L+D 


Therefore the given system of linear equations in Eqn.(l.l) 
can be rewritten as: 

(L* +U)X=b 

Using this we get 

L*X=b-UX (1.7) 

The Gauss-Seidel method is an iterative technique that solves 
the left hand side of this expression for X, using previous 
values for X on the right hand side. Using Eqn.(l .7)iteratively 
,it can be written as: 

X«= Lr^b-UX^-^), k=l,2,...., (1.8) 

provided that L. exists. Substituting for L* from above we 
get 

X« = (D + Ly'CbMD+LVuX^-V 
k=i,2,:..„ (i.9) 

where T = -((D + L) _1 )U and C=(D + L) 'b 
Therefore Gauss-Seidel technique has the form 

X^TX^-^+C, k=l,2,.... (1.10) 

However ,by taking advantage of the triangular forms 
of D,L,U and as in Jacobi's 

method ,the elements of X k can be computed sequentially 
by forward substitution using the following form of Eqn. 
(1.5) above: 

tk-ti 
>i%j Xj 

ik-ij 

E^j *3 

i=l,2,...,n , k=l,2,3,. 

The computation ofXj ' uses only the previous elements of 

X k “ : that have already been computed in advanced to 
iteration k.This means that ,unlike the Jacobi method, only 
one storage vector is required as elements can be over 
rewritten as they are computed ,which can considered as an 
advantageous for large problems .The computation for each 
element cannot be done in parallel. Furthermore, the values at 
each iteration are dependent on the order of the original 
equations. 

E. The Matlab Program for Gauss-Seidel iteration 
method 

The Matlab Program for Gauss-Seidel Method with it’s 
Command Window is shown in the Fig.(1.2) 




F. SOR iteration method (Successive Over Relaxation) 

SOR method is devised by applying an extrapolation w to 
the Gauss-Seidel method .This extrapolation takes the form of 
a weighted average between the previous iteration and the 
computed Gauss-Seidel iteration successively. 

If w>0 is a constant, the system of linear equations in 
Eqn.(1.1) can be written as 

(D+wL)X=wb-[wU+(w-1 )D] X (1.12) 

Using this Eqn.(1.12) the SOR iteration method is given by 


X :k = (D + wL) -1 [(l - w)D- wU]X tk - D 
+w(D + wL) _L b k=l,2,3... (1.13) 
provided'-.D T wZJ -1 exists . 


function jacobi2{A.b.tol.x.N) 

%jaoobi(A. b. tol.X.K) solve iteratively a system of linear equations %whereby 
%A is tli? coefficient matrix, and b is the right-hand side column vector 
%tol : relative residual error tolerance for break condition 


%X : Nxl start vector (the initial guess) 


%N is the maximum numb er of iterations 

%The method implemented is the Jacobi iterative 

%The starting vector is the null vector, but can be adjusted to one's needs 

%The iterative form is based on the Jacobi transition iteration matrix 


%Tj = inv(D)*(L+TJ] and the constant vector cj = mv(D) ^b 
%The output is the solution vector X 
%sphtting matrix A into the three matrices L, L 1 and D 
D = diag{diag(A)); 

L = tril(-A.-l); 

U= mu(-AJ); 

%transition matrix and constant vector us ed for iterations 
Tj = inv(D)*(L+U); 
cj = mv(D) ?: b: 
k=l; 

while k - r -= N 


» A=[0.2 0.11 1 0 ; 0.1 4-11 -1: 
1-1600-2; 1 1 0 S 4 : 
0-1-2 4 700]; 
»x=[0;0;0;0;0]; 

» b=[l;2^;4p]; 

» toH) 00005 ; 

» N=91 ; 

»jacobi2(A .b .tol.x.K) 

The procedure was successful 
Condition x'"(k+l) - x" (k) - r -tol was 
met after k iterations 
91 
x= 


x(: .k+ll) = Tj *x(: .k) + cj ; 

B= norm(x( ;k+ 1 )-x( :.k)); 
if B- r -tol 

dispCTheprocedure was successful') 


7.8597 

0.4229 

-0.0736 

-0.5406 

00106 


dispCCenditisn Hx J \k+l) - -"-tol was met after k iterations; 


disp(k); disp\ :-i = ');di3p(x(:.k— 1)>; 
break 


end 


k = k+l; 
end 

if B=-tol k =■ N 

diapCN'IaKiinuin number of iterations readied without satisfying conditionf);disp(C k'Xk+1) - x' Xk) 
-^toO; dispito!); 

dispCPlease. examine the sequence of iterates r ) 

dispCIn case yon obser^’e convergence, then increase the maximum number of iterations') 
disp<x); 

end _ 

rig.<l.l) The JACOBI Matlab Pro Siam with its 
Command Window for RUNNING it. 


function gauss_seidel(A ! b : toLx, N) 

%Gauss_Seidel(A b : tol : X, N) solve iteratively a system of linear 0 oequationswhereby 
%A is the coefficient matrix, and b is the right-hand side column 
%vector 

%tol: relative residual error tolerance for break condition 
%X: Nxl start vector (the initial guess) 

%N is the maximum number of iterations 
%The method implemented is the Gauss-Seidel iterative 
%The starting vector is the null vector, but can be adjusted to one's needs 
°/oThe iterative form is based on the Gauss-Seidel transition, iteration matrix 
%Tg = inv(D-L)*U and the constant vector eg = inv(D-L)*b 
%The output is the solution vector x 
%splitting matrix A into the three matrices L, U and D 


D = diag(diag(A)): 

L = tril(-A,-1); 

U = triu(-A,l); 

°/otransition matrix and constant vector used for iterations 
Tg = inv(D-L)*U; 
cg = inv(D-L)*b: 
k = 1; 

while k <= N 
x(:,k-l) = Tg*x(:,k) + eg: 

B=n orm(x(: ,k-1 )-x(: ,k)); 

I f B<tol 

disp(The procedure was successful') 
disp('Condition ||x / '<k—1) - x"{k)| <tol 
was met after k iterations') 
disp(k); disp(’x = ');disp(x(:Jcvl)); 
break 
end 

k = k—1: 
end 


» A= [0.2 0.1 1 1 0 ; 0.1 4-1 1-1:1-1 60 0-2: 
1 1 0 S 4 ; 0 -1 -2 4 700]; 


» x= [0:0:0:0:0]: 

» b= [1:2;3:4:5]: 

» N=32; 

» tol=0.00005 ; 

» gauss_seidel(A,b,tol,x,N) 

The procedure was successful 
Condition ||x'"fk-l) - x'{k)|| <tol 
was met after k iterations 31 


x = 


7.8596 

0.4229 

-0.0736 

-0.5406 

0.0106 


if B>tol || k> N 

disp('Maximum number of iterations readied without satisbing condition:') 
dispCHx'Tk-l) - x “-(k)|| <tol'); disp(tol): 
disp('Please, examine the sequence of iterates') 

disp('In case you obsewe convergence, then increase the maximum number of iterations') 

disp(x): 

end 


% &&&&& IN THE COMMAND WINDOW TYPE TO RUN THE PROGRAM &&&&&&&& 
% A=[0.2 0.1 1 10:0.14-1 1-1; 1 -1 60 0-2; 1 1 0 8 4:0-1-24700] 

% N=20 
% b=[l:2:3:4:5] 

% x=[0;0:0:0:0] 

% tol=0.00005 
% gauss_seidel(A 


Fig.(1.2) The GailSS-Seidel Matlab Program with its 
Command Window for RUNNING it. 

b,tol,x, N) 
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function S0R(A : b : tol : x : N) 

%sor(A, b : tol,x,N) solve iteratively a system of linear equations %whereby 
%A is the coefficient matrix : and b is the right-hand side column vector 
%tol: relative residual error tolerance for break condition 
?bX: Nxl start vector (the initial guess) 

%N is the maximum number of iterations 
%The method implemented is that of Successive Over Relaxation 
%The starting vector is the null vector, but can be adjusted to one's needs 
%The iterative form is based on the SOR. transition iteration matrix 


%Tw = inv(D-\v*L)*((1 -w)*D+w*U) and the constant vector 
% cw = w*inv(D-w*L)*b 

%The optimal parameter wis calculated using the spectral radius of the 
%Jacobi transition matrix Tj = inv(D)*(L+U) 

%The output is the solution vector X 


%splitting matrix A into the three matrices L, U and D 
D = diag(diag(A)); 

L = tril(-A,-l); 

U = triu(-A : l); 

Tj = inv(D)*(L+U); °aJacobi iteration matrix 
rho_Tj = max(abs(eig(Tj))); 0 ospectral radius of Tj 
w= 1.25;®ooverrelaxation parameter 
dispfw - );disp(w); 

T\v = inv(D-w*L)*((l-w)*D+w*U); %SOR iteration matrix 
cw = w*inv(D-w*L)*b; ° oconstant vector needed for iterations 

k= 1; 

while k <= N 

x(:,k+l) = Tw*x(:,k) + cw; 

B=norm(x(:,k+l)-x(:,k)); 
if B<tol 

disp(The procedure was successful') 

dispCCondition ||x A (k+l) - x A (k)|| <tol was met after k iterations') 
disp(k); dispfx = ');disp(x(:,k+l)); 
break 
end 

k = k+l; 
end 

ifB>tol || k>=N 

dispf Maximum number of iterations reached without satisfying condition:') 
dispC||x A (k+l)- x A (k)|| <tol'); disp(tol); 
dispCPlease : examine the sequence of iterates') 

dispfln case you observe convergence, then increase the maximum number of iterations'); disp(x); 
end 

% &&&&& IN THE COMMAND WINDOW TYPE TO RUN THE PROGRAM &&&&&& 

% A=[0-2 0.1 1 1 0 ; 0.1 4 -1 1-1; 1 -1 60 0 -2; 1 1 0 8 4; 0-1 -2 4 700]; 

% N=30; 

% b=[U;3;4;5]; % x=[0;0;0;0;0]; 


»A=[0.2 0.1 1 1 0 ; 0.1 4 -1 1-1; 1-1 60 
0-2; 1 1 0 8 4; 0 -1 -2 4 700]; 

»N=32; 

»tol=0.00005; 

»x=[0;0;0;0;0]; 

»b=[l‘2;3;4;5]; 

» SOR(A, b : tol,x, N) 
w= 1.2500 

The procedure was successful 
Condition ||x A (k+l)- x A (k)|| <tol 
was met after k iterations 15 
x = 

7.8597 

0.4229 

-0.0736 

-0.5406 

0.0106 


Fig.(1.3) The SOR Matlab Program with its 

Command Window for RUNNING it. 


% tol=0.00005 ; % SOR(A, b,tol,x, N) 


G. The Matlab Program for SOR iteration method 

The Matlab Program for SOR Method with it it’s 
Command Window is shown in Fig.(1.3) 


II. The conjugate gradient method 

It is used to solve the system of linear equations 

Ax = b (2.1) 

for the vector x where the known n-by-n matrix A is 
symmetric (i.e. A T = A), positive definite (i.e. x T Ax > 0 for all 
non-zero vectors x in R n ), and real, and b is known as well. 
We denote the unique solution of this system by x*. 

A. The conjugate gradient method as a direct method 

We say that two non-zero vectors u and v are 

conjugate (with respect to A) if 

T 

11 Av = 0. (2.2) 

Since A is symmetric and positive definite, the left-hand side 
defines an inner product 

(U,V)a : = (AU,V) = (U,A l V) = (U.. AV) = U 1 AV (2.3) 

Two vectors are conjugate if they are orthogonal with respect 
to this inner product. Being conjugate is a symmetric relation: 
if u is conjugate to v, then v is conjugate to u. (Note: This 
notion of conjugate is not related to the notion of complex 
conjugate.) 

Suppose that {} is a sequence of n mutually conjugate 


Now , let T w = (D +- wL) -L [(l - wj D - wU] 
and C w = wtD + wL] -1 b 
Then the SOR technique has the iterative form 

X r " - T w X ! " k_1 - -f C^ r for a constant w> 0, 
k=l,2,..., (1.14) 

By Eqn. (1.14) it is obvious that the method of SOR 


is an iterative technique that solves the left hand side of this 
expression for X ^ when the previous values for X k ”~ on the 
right hand side are computed .Analytically ,this may be 
written as: 

X ;:ir - = (D + wL) _1 

(wb - [wU + Cw — ljD]X k-1 -), 
k=l,2,3... (1.15) 

By taking advantage of the triangular form of 

(D+ wL) it can proved that CD + wO "" = D 

Using this, the elements of X v can be computed sequentially 

using forward substitution as in Eqn.(1.3) in section l.land 

Eqn.(l.lO) in section 1.2 above , i.e. 

Jk) f., \ tk— l) 

Xj = U - wjxj 

-(-^(bj - xf -1 - Sjcjaij x, k_1 ), 

i=l,2,...,n, k=0,l,.... (1-16) 

The choice of the relaxation factor w is not necessarily easy 
and depends upon the properties of the coefficient matrix .For 
positive -definite matrices it can be proved that 0<w<2will 
lead to convergence, but we are generally interested in faster 
convergence rather than just convergence. 


directions. Then the{ } form a basis of R”, so we can expand 
the solution x^of Ax = b in this basis: 

n 

X. = y Kj Pj (2.4) 

i =1 

and we see that 

n 

b = AX. = y Kj A ¥> (2.5) 

i =i 

The coefficients are given by 

n 

$b = f$AX. = T Kj Fg AF] = (2.6) 

1 = 1 


(because Vi ^ AX P . : , P- are mutually conjugate) 



P^A P K 


<P'K> bi 

(PR - PK>4 


(PK ^ b) 
IIPlCllJ 


(2.7) 


This result is perhaps most transparent by considering the 
inner product defined above. 

This gives the following method for solving the equation Ax = 
b: find a sequence of n conjugate directions, and then 
compute the coefficients °-k. 

B. The conjugate gradient method as an iterative method 

The direct algorithm was then modified to obtain an 
algorithm which only requires storage of the last two residue 
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vectors and the last search direction, and only one matrix 
vector multiplication. 

The algorithm is detailed below for solving Ax = b 
where A is a real, symmetric, positive-definite matrix. The 
input vector x 0 can be an approximate initial solution or 0. 

:= b - Ax : 

P E i 

K := 0 


repeat 


P^AP k 


*Jf+l ; = p k 

r K + l := r K — Pk 

if r x +1 is sufficiently small then exit loop 

T 

^k+i^+i 


Pk : = 

fy+i ;= Fk+l + 
k := k + 1 


T 

f’K 


end repeat 
The result is *#+1 


This is the most commonly used algorithm. The same formula 
for pn is also used in the Fletcher-Reeves nonlinear 
conjugate gradient method. 


C. The code in Matlab for the Conjugate 

Iterative Algorithm in Fiq.(1.4) & Fiq.(1.5) 

Alternatively another Matlab Code for Conjugate 
Gradient Algorithm is in Fiq.(1.5) 


D. Convergence properties of the conjugate gradient 
method 

The conjugate gradient method can theoretically be viewed 
as a direct method, as it produces the exact solution after a 
finite number of iterations, which is not larger than the size of 
the matrix, in the absence of round-off error. However, the 
conjugate gradient method is unstable with respect to even 
small perturbations, e.g., most directions are not in practice 
conjugate, and the exact solution is never obtained. 
Fortunately, the conjugate gradient method can be used as an 
iterative method as it provides monotonically improving 
approximations x- K to the exact solution, which may reach 

the required tolerance after a relatively small (compared to the 
problem size) number of iterations. The improvement is 
typically linear and its speed is determined by the condition 
number K (A j of the system matrix A: the larger is K (A'i, the 

slower the improvement. 

If K(A) is large, preconditioning is used to 


E. The preconditioned conjugate gradient method 

Preconditioning is necessary to ensure fast convergence 
of the conjugate gradient method. The preconditioned 
conjugate gradient method takes the following form: 

r : := b - Axfl 
z c , := 

P E , := 

K := 0 




repeat 

r 


P l-A P k 


+ i X K Pk 

r ff+l := r K Pk 

if r k+ i is sufficiently small then exit loop end if 
Zk + l : = M “ lr k+1 


Pk '= 


T 

s k+i r k+i 

T 

z k rk 


P« + L ; = Zk + 1 + P k P k 

k :=k+l 


end repeat 

The result is 

The above formulation is equivalent to applying the 
conjugate gradient method without preconditioning to the 
system 

E - 1 A(E -1 -) t x = E _1 b 

where EE" = M and x= E"k 


The preconditioned matrix M has to be symmetric 
positive-definite and fixed, i.e., cannot change from iteration 
to iteration. If any of these assumptions on the preconditioned 
is violated, the behavior of the preconditioned conjugate 
gradient method may become unpredictable. 

An example of a commonly used preconditioned is the 
incomplete Cholesky factorization. 


F. The code in Matlab for the Preconditioned Conjugate 
Gradient Algorithm in Fig(1.6) 


G. The flexible preconditioned conjugate gradient 
method 


In numerically challenging applications 
sophisticated preconditioners are used, which may 
lead to variable preconditioning, changing between 
iterations. Even if the preconditioned matrix is 
symmetric positive-definite on every iteration, the 
fact that it may change makes the arguments above invalid, 
and in practical tests leads to a significant slowdown of the 
convergence of the algorithm presented above. Using the 


Polak-Ribiere formula 



r K 


replace the original system 

Ax - b = 0 with M -1 (Ax - b) = 0 so that 

K gets smaller than K(A) , see below 


instead of the Fletcher-Reeves formula 


Pk 




may dramatically improve the convergence in this case. This 

version of the preconditioned conjugate 

gradient method can be called flexible, as it allows 
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function [x, niter, flag] = CONJUGATE_GRAD(A, b, s, tol, maxiter) 
°/oCG Conjugate Gradients method 
°/bInput parameters 

°/oA : Symmetric, positive definite NxN matrix 
°/ob : Right-hand side Nxl column vector 
°/os : Nxl start vector (the initial guess) 

°/otol: relative residual error tolerance for break condition 
°/omaxiter : Maximum number of iterations to perform 


% Output parameters 

%x : Nxl solution vector 

%niter : Number of iterations performed 

%flag: 1 if convergence criteria specified by TOL could 

%notbe fulfilled within (the specified maximum Successful) 

%number of iterations, 0 otherwise 
x= s; % Set xO to the start vector s 
r = b - A*s; % Compute first residuum 

p = r; 

rho = r'*r; 

niter = 0; % Init counter for number of iterations 
flag = 0; % Init break flag 

%Compute norm of right-hand side to take relative residuum as 
%break condition 
normb = norm(b); 

if normb<eps % if the norm is very close to zero, take the 
%absolute residuum instead as break condition 
%norm(r) >tol since the relative 

warning(['norm(b) is very close to zero, taking absolute as a break condition']); 

normb = 1; 

end 

while (norm(r)/normb>toI) % Test break condition 
a = A*p; 


alpha = rho/(a'*p); 
x = x + alpha*p; 
r = r - alpha^a; 
rhonew = r'*r; 
p = r + rho_new/rho * p; 
rho = rhonew; 
niter = niter + 1; 

% if max. number of iterations 

if (niter = maxiter) 

flag = 1; % is reached, break 

break 

end 

end 


» A=[0.2 0.1 11 0 : 0.1 4 -11 -1; 1 -1 60 0 -2; 1 1 0 8 4 ; 0 -1 -2 4 700]; 
» maxiter=6; 

»b=[l;2;3;4;5]; 

» s=[0;0;0;0;0]; 

»tol=0.00005; 

» [x, niter, flag] = CONJUGATE_GRADIENT(A, b, s, tol, maxiter) 


x = 7.8597 0.4229 
niter = 5 
flag = 0 


-0.0736 -0.5406 0.0106 


% &&&&& IN THE COMMAND WINDOW TYPE TO RUN THE PROGRAM &&&&&& 
% A=[0.2 0.1 1 1 0 ; 0.1 4 -1 1 -1; 1 -1 60 0 -2; 1 1 0 8 4 ; 0 -1 -2 4 700]; 

% maxiter=6; - 

% b=[l;2;3;4;5]; Fig.(1.4) The CONJUGATE_GRADIENT Matlab Program with its 

% s=[0;0;0;0;0]; Command Window for RUNNING it & below it the output 

% tol=0.00005; 

% [x, niter, flag] = CONJUGATE_GRADIENT(A, b, s, toL maxiter) 


particular, it converges not slower than the locally optimal 
steepest descent method. 


III. Conclusion 

When We ran each Matlab Code program and testified 
with a numerical example We observe that 
the non-stationary algorithms converge faster than the 
stationary algorithms. 

The reader can look and compare at the output result 
programs in Figures denoted above ((these Figures are 
Fig.(l) , Fig(2) , ....,Fig(5) and Fig(6) )) to ensure that 
non-stationary algorithms especially 

the Preconditioned Conjugate Algorithm approximates the 
solution accurately to five decimal 

places with the number of iterations N equals 5 which less 
compared to the other Algorithms tackled in this 
paper. 

The paper also analysis the Convergence properties of the 
conjugate gradient method and shows there are two methods, 
the direct or the iterative conjugate. 

Further, better convergence behavior can be achieved 
when using Polak-Ribiere formula which 
converges locally optimal not slower than the locally 
optimal steepest descent methed. 

Furthermore, these Algorithms can are 
recommended for similar situations which are arise in many 
important settings such as finite differences, 
finite element methods for solving partial differential 
equations and structural and circuit analysis. 


function [iterationno , x] = conjgrad{A : b,x) 
r=b-A*x; 


p=r: 

rsold=r'*r: 


»A=fl>-2Q-l 110:0.1 4 -1 1-1: 1-1600-2: 1 10 E4 : 0-1-2 
4 700]; 


for i=l:10000000 
Ap=A*p; 

atpha=r sold ''(p 1 * Ap) ; 
x=x—alpha*p: 
r=r-aIpha*Ap; 


»b=[l;2;3;4;5]; 

» x=[0;0;0;0;0]; 

» [iterationno , x] = conjgrad{A ; b : x) 


iterationno =6 


rsnew=r'*r; 
if sqrt(rsnew)<5e-10 
break: 


x= 7.8597 
0.4229 


end 

p=r+rsnew r sold*p: 
rsold=rsnew; 


-0.073d 

-0.5406 

0.0106 


end 

iterationno=i 


Fig. (1.5) T fa e CONJUGATEGRADIENT Ma da b Pro gi am 
its Command Window for RUNNING it 
& below if the output result 


% & IN THE COMMAND WINDOW TYPE TO RUN THE PROGRAM 


with 


& 


% A=[0.2 0.1 110: 0.1 4 -1 1 -1; 1 -1 60 0 -2; 1 1 0 8 4 ; 0 -1 -2 4 700]: 
% b=[l;2;3;4;5]; 

% x=[0;0;0:0;0]: 

% [iterationno , x] = conjgrad{A : b ; x) 


function Preconditioned_Conjugate(A,b : c : toLN) 
[n.m]=size(A); 
invc=inv(c); 
x=zeros(n.l); 
r=b; 

p=invc*b; 
y=mvc*r; 
eir(l)=norm(r); 
k=0; 

while nonn(r) >tol<fc&k<N 
k=k+l; 
z=A*p; 

alpha=(y'*r)' , (p’*z); 
x=x-alpha*p: 
mew=r-alpha *z; 
yne w=inv c * mew; 
beta=(ynew'*mew')(y’*r); 
p=ynew-b eta*p; 
r=mew; 


» A=(0.2 0.111 0 ; 0.1 4 -11 -1; 1 -160 0 -2; 1 1 0 8 4 ; 0 -1 -2 4 700] 

>>N=5 ; 

» b=[l;2;3;4;5]; 

»c=eye(5,5); 

»tol=0.00005 ; 

» Preconditioned_Conjugate(A : b,c.toLN) 

x= 7.8597 0.4229 -0.0736 -0.5406 0.0106 
err= 7.5271 5.5600 0.7239 0.5572 0.0000 


Fig.(1.6) The Preconditioned Conjugate Matlab Program with 
Command Window for RUNNING it & below it the output result 


v=vnew; 
err(k)=nonn(r); 
end 

dispC x=0; disp(x'); 
dispf en=’); disp(eir); 

% &&&&& IN THE COMMAND WINDOW TYPE TO RUN THE PROGRAM 
&&&&&&&& 


% A=[0.2 0.1 1 1 0 ; 0.1 4 -1 1 -1; 1 -1 60 0 -2; 1 1 0 8 4 ; 0 -1 -2 4 700]; 
% N=5 ; 

% b=[l;2;3;4;5]; 

% c=eve(5,5); 

% tol=0.00005; 

% Preconditioned_Conjugate(Ab.c.toLN) 


for variable preconditioning. The implementation of 
the flexible version requires storing an extra vector. 

T 

For a fixed preconditioned, ^+^=0 so both 

formulas for are equivalent in exact arithmetic, i.e., 

without the round-off error. 

The mathematical explanation of the better 
convergence behavior of the method with the Polak-Ribiere 
formula is that the method is locally optimal in this case, in 
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